import matplotlib.pyplot as plt
from scipy.stats import scoreatpercentile
import numpy as np
import csv
import random
from datetime import datetime
from time import strftime
from sympy import *

'''
revision of test24
to calculate the false nagetive, false positive

'''

def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV):
    sKey = []
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        for item in ra.fieldnames:
            temp = float(record[item])
            sKey.append(temp)
    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def fivenum(v):
    """Returns Tukey's five number summary (minimum, lower-hinge, median, upper-hinge, maximum) for the input vector, a list or array of numbers based on 1.5 times the interquartile distance"""
    import numpy as np
    from scipy.stats import scoreatpercentile
    try:
        np.sum(v)
    except TypeError:
        print('Error: you must provide a list or array of only numbers')
    q1 = scoreatpercentile(v,25)
    q3 = scoreatpercentile(v,75)
    md = np.median(v)
    return np.min(v), q1, md, q3, np.max(v),

def nullH():
    c = 6000
    #clusterPop = [1684327, 360275, 7627173]
    clusterPop = [1611198, 501040, 7025156]
    
    #n = clusterPop[2]
    N = 29535210
    #k = 1.645 * ((c*n*(N-n)/(N*N))**0.5) + c*n/N
    print '#CASES = ', c
    print 'average risk = ', (c + 0.0)/totalPop
    risk = [(c + 0.0)/totalPop]
    for n in clusterPop:
        print '----------'
        print 'pop = ', n
        k = Symbol('k')
        k = solve((k-c*n/N)/((c*n*(N-n)/(N*N))**0.5) - 1.645, k)
        k = k[0]
        r = Symbol('r')
        eqn = Eq((((N - n + n * r) * k - c * n * r) ** 2)/(c * n * r * (N - n)), 3.09 ** 2)
        r = 0
        temp = solve(eqn)
        for i in temp:
            if i > 1:
                r = i
        if r > 0:
            print 'r = ', r
        else:
            print 'error in solving r'

        print 'E(c|Ha) = ', c * n * r / (N - n + n * r)
        print 'Var(c|Ha) = ', c * n * r * (N - n) / ((N - n + n * r) ** 2)
        print 'Var(c|Ha) = ', c * n * r * (N - n + n * r - c * r) / ((N - n + n * r) ** 2)
        print 'risk = ', c * r / (N - n + n * r)
        risk.append(c * r / (N - n + n * r))
    print risk

#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print "begin at " + getCurTime()
    #print np.__version__
    # total pop: 29535210

    '''
    # old risk area
    mixed = [91,98,101,104,114,115,119,126,131,142,146,147,154,162,168,172] # pop: 1684327
    rural = [8,9,10,11,12,13,14,15,17,19,20,26,28,33,34,37] # pop: 360275
    urban = [105,107,112,120,122,125,127,128,130,133,134,141,143,149,152,155]   # pop: 7627173
    '''
    # new 2 risk area
    mixed = [95,99,109,110,114,115,117,119,126,131,139,140,142,146,147,237] # pop: 1611198
    rural = [9,10,12,13,14,17,20,23,26,27,28,33,34,35,38,41]    # pop: 501040  
    urban = [125,127,128,130,133,134,136,138,141,143,149,151,152,156,157,163] # pop: 7025156

    hot_1 = [9,130,147]
    hot_2 = hot_1 + [10,133,154]
    hot_4 = hot_2 + [12,17,125,131,141,146]
    hot_8 = hot_4 + [14,19,20,26,114,115,119,120,128,134,149,168]
    hot_16 = mixed + rural + urban
    
    #unitCSV = "C:/TP1000_1m.csv"
    #unitCSV = 'C:/_DATA/CancerData/SatScan/mult6000/three16_format.csv'
    unitCSV = 'C:/_DATA/CancerData/SatScan/own/6000.csv'
    dataMatrix = build_data_list(unitCSV)  # [id, pop, cancer1, cancer2, cancer3]
    unit_attri = np.zeros((len(dataMatrix),4))
    
    # unit_attri [cancer, pop, rate, id]
    unit_attri[:,1] = dataMatrix[:,1]
    unit_attri[:,0] = dataMatrix[:,2]

    totalPop = 0
    for item in unit_attri[:,1]:
        totalPop += item

    pop = 0
    for item in urban:
        pop += unit_attri[item,1]

    cancerT = []
    for c in range(0,1000):
        cancer = 0
        for item in mixed:
            cancer += dataMatrix[item,2 + c]
        cancerT.append(cancer)
    #print np.var(cancerT)
    #print np.mean(cancerT)
    print fivenum(cancerT)
        
    #print pop*6000*1.15/(totalPop-pop+pop*1.15)

    #nullH()
    
    print "end at " + getCurTime()
    print "========================================================================"  

           
